###################
###################
###Matching Code###
###################
###################
library(MatchIt)
library(marginaleffects)
library(Hmisc) 
library(plotrix) 
set.seed(1234)
###Match on key issues
full.1996.2019$outbreak_bin <- ifelse(full.1996.2019$outbreak>0,1,0)
##Remove NA data
subset.dat <- na.omit(full.1996.2019[,c("acled_battle_state","acled_battle_rebel","acled_battle_milall","acled_battle_polmil","acled_battle_idmil","log.nl","logged.pop","lagacled_battle_state","lagacled_battle_rebel","lagacled_battle_milall","lagacled_battle_polmil","lagacled_battle_idmil","t","spei","p_anom","outbreak_bin","t_anom","gid","month", "year")])  
#Create a logical binary 
subset.dat$f_outbreak_bin <- as.logical(subset.dat$outbreak_bin)

####################
###State Conflict###
####################
#Run PSM algorithm
match.dat <- matchit(f_outbreak_bin ~  log.nl + logged.pop +
                       lagacled_battle_state +
                       p_anom+spei+t_anom+month, 
                     data = subset.dat, method="nearest", ratio=1)
#summary(match.dat)
#Turn into a dataframe
match.dat2 <- match.data(match.dat)
#Fit linear model
fit1 <- lm(acled_battle_state ~ 
             f_outbreak_bin*(log.nl + logged.pop +
                               lagacled_battle_state + 
                               p_anom+spei+t_anom),
           data = match.dat2, weights = weights)
summary(fit1)
#Estimate ATT for outbreak binary
marginaleffects::avg_comparisons(fit1, variables = "f_outbreak_bin",
                                 newdata = subset(match.dat2, f_outbreak_bin == 1),
                                 wts="weights")

####################
###Rebel Conflict###
####################
#Run PSM algorithm
match.dat21 <- matchit(f_outbreak_bin ~  log.nl + logged.pop +
                         lagacled_battle_rebel +
                         p_anom+spei+t_anom+month, 
                       data = subset.dat, method="nearest", ratio=1)
#summary(match.dat21)
#Turn into a dataframe
match.dat22 <- match.data(match.dat21)
#Fit linear model
fit2 <- lm(acled_battle_rebel ~ 
             f_outbreak_bin*(log.nl + logged.pop +
                               lagacled_battle_rebel + 
                               p_anom+spei+t_anom),
           data = match.dat22, weights = weights)
summary(fit2)
#Estimate ATT for outbreak binary
marginaleffects::avg_comparisons(fit2, variables = "f_outbreak_bin",
                                 newdata = subset(match.dat22, f_outbreak_bin == 1),
                                 wts="weights")

#################################
###Political Militias Conflict###
#################################
#Run PSM algorithm
match.dat31 <- matchit(f_outbreak_bin ~  log.nl + logged.pop +
                         lagacled_battle_polmil +
                         p_anom+spei+t_anom+month, 
                       data = subset.dat, method="nearest", ratio=1)
#summary(match.dat21)
#Turn into a dataframe
match.dat32 <- match.data(match.dat31)
#Fit linear model
fit3 <- lm(acled_battle_rebel ~ 
             f_outbreak_bin*(log.nl + logged.pop +
                               lagacled_battle_polmil + 
                               p_anom+spei+t_anom),
           data = match.dat32, weights = weights)
summary(fit3)
#Estimate ATT for outbreak binary
marginaleffects::avg_comparisons(fit3, variables = "f_outbreak_bin",
                                 newdata = subset(match.dat32, f_outbreak_bin == 1),
                                 wts="weights")

#################################
###Identity  Militias Conflict###
#################################
#Run PSM algorithm
###Match on key issues
match.dat4 <- matchit(f_outbreak_bin ~  log.nl + logged.pop +
                        lagacled_battle_idmil +
                        p_anom+spei+t_anom+month, 
                      data = subset.dat, method="nearest", ratio=1)
#summary(match.dat4)
#Turn into a dataframe
match.dat42 <- match.data(match.dat4)
#Fit linear model
fit4 <- lm(acled_battle_idmil ~ 
             f_outbreak_bin*(log.nl + logged.pop +
                               lagacled_battle_idmil + 
                               p_anom+spei+t_anom),
           data = match.dat42, weights = weights)
summary(fit4)
#Estimate ATT for outbreak binary
marginaleffects::avg_comparisons(fit4, variables = "f_outbreak_bin",
                                 newdata = subset(match.dat42, f_outbreak_bin == 1),
                                 wts="weights")

######################
###Plot the Results###
######################

##Barplots for each group
#The means and 95% CI for each categroy
state.change <- c((-0.0004), (-0.021), 0.020)
reb.change <- c(0.005, (-0.008), 0.017)
milpol.change <- c(0.0002, (-0.017), 0.018)
milid.change <- c(0.020, 0.003, 0.036)
change.dat <- data.frame(rbind(state.change,reb.change,milpol.change,milid.change))
names(change.dat) <- c("mean", "lci", "hci")
#Plot the bars
pdf("Fig2.pdf", width = 8, height = 6)
df.bar <- barplot(change.dat$mean, col=c("darkblue","orange","purple","yellow"), ylab="ATT",
                  width=0.35, xlim=c(0,1.75), ylim=c((-0.03),0.055),
                  names.arg=c("State", "Rebels", "Pol. mil.", "Id. mil.")) 
errbar(df.bar[,1], change.dat$mean, yplus=change.dat$hci, yminus =change.dat$lci, add=T, xlab="")
#Add legend
legend(x = "topright",          # Position
       legend = c("Bar plots represent mean values with +/- 1.97 SEM"))   # Legend texts

dev.off()

